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Summary. We report a three-variable simplified model of excitation fronts in hu- 
man atrial tissue. The model is derived by novel asymptotic techniques from the 
biophysically realistic model of Courtemanche et al. [11] in extension of our previ- 
ous similar models. An iterative analytical solution of the model is presented which 
is in excellent quantitative agreement with the realistic model. It opens new possi- 
bilities for analytical studies as well as for efficient numerical simulation of this and 
other cardiac models of similar structure. 
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1.1 Introduction 

The mechanical activity of the heart is controlled by electrical impulses prop- 
agating regularly through it during our entire lifespans [22] . A disturbance in 
the regular propagation may lead to life-threatening cardiac arrhythmias [31]. 
Sudden cardiac death, for instance, accounts for 300 000 to 400 000 deaths 
annually in the United States alone [14, 23], i.e. more than AIDS, breast 
and lung cancer. This entails intensive research into the mechanisms of heart 
functioning and failure. The accumulated information reveals an overwhelm- 
ing complexity of the patterns of electrical cardiac activity. 

True understanding of the experimental data requires the development of 
appropriate qualitative cardiac models [10, 19, 20]. One approach to cardiac 
modelling is to take into account the various levels of membrane, cellular and 
myocardial structure and their interactions and to model the action potential 
(AP) on the basis of experimental measurements of ion fluxes in as much de- 
tail as possible. The resulting models are known as realistic or detailed ionic 
models. The first example of this type of models was developed by Noble 
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[25, 26] and there are now such models for various cardiac cells in different 
species, e.g. [3, 12, 21, 27, 35] and many others. However, since these models 
are very complex and highly nonlinear, it is difficult to assess contribution 
of specific model components to different patterns of activity. Furthermore, 
their computation is arduous because they contain a large number of equa- 
tions and small parameters. They become very expensive and time-consuming 
especially when large volumes of tissue are simulated. One possible alternative 
is to search for simplified models which could mimic the most important AP 
properties, allow analytical studies and reduce the computing requirements. 
Many simplified models have been suggested, either phenomenologically, or 
based on the structure of the realistic models, [1, 4, 13, 15, 18] etc.. However, 
all of these models contain arbitrary elements in the sense that they are not 
derived from any of the realistic biophysical cell models and lack explicit 
correspondence with the biophysical structure of the cardiac tissues. E.g. 
van der Pol and van der Mark modelled heartbeat as an electronic relaxation 
oscillator [28]. 

Our recent studies [8, 32, 33, 7] have demonstrated a serious disadvan- 
tage of the most popular and successful simplified generic model of cardiac 
excitability ~ the FitzHugh-Nagumo (FHN) equations [17, 24]. It cannot de- 
scribe adequately the one feature of excitation propagation which is most 
important for medical applications, namely the way regular propagation fails 
and arrhythmias occur. An excitation wave in real atrial tissue, or in a realistic 
model, may fail to propagate if the temporal gradient of the transmembrane 
voltage at the front becomes too small to excite the tissue ahead of it e.g. 
if the wave fails to propagate fast enough [8]. Then the wave front looses 
its sharp spatial gradient and its further spread is purely diffusive, i.e. the 
front dissipates. This happens long before the back of the excitation impulse 
catches up with the front. This type of propagation failure does not exist in 
a FitzHugh-Nagumo type model [6] since it is known that the propagation of 
a wave front in this system may be slowed down, halted or even reversed [16]. 
This phenomenon is illustrated in figure 1.1: a temporary block of excitability 
only temporary halts the FHN wave, but completely disrupts propagation in 
the realistic model, even though it lasts much shorter than the AP. 

Earlier we have proposed novel asymptotic methods of reduction of car- 
diac equations [9]. In this paper we use these methods to derive a three- 
component description of the propagating excitation fronts and their dissi- 
pation. The virtue of our model is that it reproduces propagation failure 
unlike the FitzHugh-Nagumo models because it is derived in a reliable way 
from the realistic ionic model. We also report an analytical solution for this 
three-component model. The analytical solution is constructed as an iterative 
procedure and it may be seen as a generalisation in which the caricature solu- 
tions presented in our earlier papers [5, 6, 30] appear as first approximations. 
This, is to our knowledge, the first analytical solution, albeit in quadratures, 
which gives a numerically accurate prediction of the front propagation velocity 
(within 16%) and its profile (within 0.7 mV) in human atrial tissue. 



Page : 2 



job: asm07 



macro : svmult . els 



date/time: 25-May-2009/16 : 06 




-40 
-80 


-40 
-80 
M 'J 
1-40 
^ -80 


-40 
-80 


-40 
-80 



t — 






t=40 




t=80 


1 ^ ^ ^ ^ 


^ 1 

t=120 


1 ^ ^ ^ ^ 


t=160 













1 Model of Atrial Excitability 
N Eq. (1.3) 



t=0 




200 400 600 



«=32 




-40 
-80 


-40 
-80 


-40 
-80 


-40 
-80 


-40 
-80 



t=0 



t=76 



t=152 



30 60 

X 



100 200 300 



Fig. 1.1. Propagation of excitation in the models of Courtemanche et al. [11] (first 
column), FitzHugh-Nagumo (second column) and in equations (1.3) (third column), 
through a temporary block of excitability, introduced by artificially reducing the 
value of a parameter representing the main excitatory ionic current responsible for 
the initiation of the front. E.g. in the three-variable model (1.3) this parameter is 
j which was decreased from the normal value of 0.9775 to 0.28 during the block. 
In FHN, propagation resumes after the block is removed; in CRN and (1.3) it does 
not. 



1.2 Mathematical formulation of the problem 

Atrial tissue model. In our study atrial tissue is a one-dimensional, ho- 
mogeneous and isotropic medium satisfying a system of reaction-diffusion 
equations 

drvi^ b ■ d\u + F{u) (1.1) 

where F{\i) is a vector defined according to the atrial single-cell realistic CRN 
model [11], u = {V,m,h, j, . . <^M?^ is the vector of all dynamic variables 
of the model and D = diag(£', 0, 0, . . . ) is the tensor of diffusion in which only 
the coefficient of the voltage V is nonzero. This simplified description focuses 
on the excitation and propagation of impulses, while ignoring the effects due 
to geometry, anisotropy and heterogeneity of a real atrium. 
Asymptotic reduction. In order to reduce the dimension and complexity 
of the problem, we perform a formal analysis of the time scales of dynamic 
variables. For the system (1.1) we define characteristic time scale functions, 

Ti{ui, . . . ,U2i) = [dFi/dui) ^ and compare their magnitudes obtained nu- 
merically for a space-clamped version of the system as shown in figure 1.2. 
The variables, whose time scales Ti are relatively small, are fast variables since 
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Fig. 1.2. Asymptotic properties of tlie atrial CRN model [11]. (a) Time scale 
functions of dynamical variables vs. time, (b) Quasistationary values of the gat- 
ing variables m and h. (c) Transmembrane voltage V as a function of time, (d) 
Main ionic currents: iNa is the fast sodium current (shown scaled by a factor 0.1), 
/in = /b,Na + /NaK + /ca,L + /b,Ca + /NaCa is the sum of all othcr luward currents and 
/out — /p,Ca + /ki + /to + /Kur + /kf + /ks + /b,K is the sum of all outward currents; 
the individual currents are described in [11]. The results are obtained for a space- 
clamped version of the model at values of the parameters as given in [11]. A typical 
AP was triggered by initialising the transmembrane voltage to non- equilibrium value 
of V = -20 mV. 

they change significantly during the upstroke of a typical AP, while all other 
variables, whose time scales arc relatively large, arc slow variables because 
they change only slightly during this period. Figure 1.2(a) demonstrates that 
the variables V , m, h, Ua, w, Oa, d are fast variables comparable with the time 
scale of the AP upstroke. 

A specific feature of system (1.1) is that of the various ionic currents in 
the system only the sodium current /Na is significantly large during the AP 
upstroke, whereas other currents are small at this stage as can be seen in 
figure 1.2(d). Secondly, the fast sodium current /Na is only large during the 
AP upstroke, and almost vanishes otherwise, because either gate m or gate 
h or both are nearly closed outside the upstroke since their quasistationary 
values rri{V) and h{V) are small there as illustrated in figure 1.2(b). 
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To formalise the distinction between fast and slow terms we perform an 
asymptotic embedding of system (1.1). We introduce an artificial parameter e 
into the system so that for e = 1 the original system is recovered, while in the 
limit e — > only the terms comparable with the time scale of the AP upstroke 
arc retained, 



drV = DdlV - 



{e-'lMV,m,h,j) + S'j{V,...)) 
'TV ^ uuxv , 

(m(y;e)-m) rm(y), e = 1, 

{y{V) - y) 

"Ty = 777^' y = Ua,w,Oa,d, 

drV ^WiV,...), (1.2) 

where 9{) is the Hcavisidc function, is the sum of all currents except 

the fast sodium current /Na, the dynamic variables V, to, h, Ua, Oa and d 
arc defined in [11], U = (j, Oi, . . . , Nai, Ki, . . . )^ is the vector of all other, 
slower variables, and W is the vector of the corresponding right-hand sides. 
Novel features of the asymptotic embedding (1.2), non-standard in comparison 
with the theory of fast-slow systems [34, 29, 2] are (a) the introduction of 
the asymptotic factor e^^ only at one term /Na ui the right-hand side of 
the equation for V whereas the standard factor e at the derivative would 
be equivalent to factor at the whole right-hand side, and (b) that in 
the limit e — > 0, functions rn(V) and h{V) have to be considered zero in 
certain overlapping intervals V G (— oo, Vm] and V £ [Vh, +oo), and Vh < Vm, 
hence the representations rn{V;Q) = 6{V — Vm) and h{V;0) = 9{Vh — V). 
These aspects, as applied to the fast sodium current, have been shown to 
be crucial for the correct description of the propagation block [5]. A more 
detailed discussion of the parameterisation (1.2) can be found in reference [9]. 

The exact value of D is not essential for the theoretical analysis, as its 
change is equivalent to rescaling of the spatial coordinate. To operate with 
dimensional velocity, we assume the values D = 0.03125, as in [8, 9] and 
Cm = Ia^F cm~^. We perform the scaling t = e^^T, x = {eD)~^/^X , take 
the limit e ^ and notice that the equations for the variables denoted by y 
in (1.2) decouple from the voltage equation. Thus we arrive at the conclusion 
that only the following three-variable system needs to be considered for a 
description of the propagation of an AP front or its failure, 

dtV ^dlV + I^iV)jhm^, (1.3a) 
dth = {e{Vh -V)- h)/Th{V), (1.3b) 
dtm = {e{V - Kn) - m)/TmiV). (1.3c) 
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In other words, we consider the fast time scale on which the upstroke of the 
AP occurs, neglect the variations of slow variables during this period as well 
as all transmembrane currents except /Na, as they do not make significant 
contribution during this period and replace m and h with zero when they are 
small. The parameters and functions in (1.3) arc defined as in [11], namely 

T^.{V)^gNa{VNa~V), (1.4a) 

Tk{V) ^ {ak{V) + pk{V))~\ k^h,m, (1.4b) 

k{V)=ak{V)/{ak{V)+pk{V)), k = h,m, 

ahiV) = 0.135 e-(^+8")/6-8 g^^y - 40), 

f3^{V) = (3.56 e" °^9^ + 3.1 x 10^ e^-^sv) g^.y ~ 40) 

+ e{V + 40) (0.13(1 + e-(^+io-66)/ii.i))-\ 

0.32(1^ + 47.13) 

am[V) - _^ _ g-o.i(y+47.i3) ' 

f3r,^{V) = 0.08e-^/", 

5jva = 7.8, VNa = 67.53, = -66.66, y„ = -32.7. 

Two new 'gate threshold' parameters 14 and Vm appear in the system and 
are chosen from the conditions h[Vh) = 1/2 and m^{Vm) = 1/2. As follows 
from the derivation, variable j, the slow inactivation gate of the fast sodium 
current, acts as a parameter of the model. It is the only one of all slow variables 
included in the vector U that affects our fast subsystem. We say that it 
describes the excitability of the tissue. 

Travelling waves and reduction to ODE. We look for solutions in the 
form of a front propagating with a constant speed and shape. So we use the 
ansatz F{z) — F{x + ct) for F = V,h,m where c is the dimensionless wave 



speed of the front, related to the dimensional speed C by c = (e/Z?)^/^C. 
Then equations (1.3) reduce to 

V" = cV' -~l^,{V)jhm^, (1.5a) 

h' ^ {cTh{V)y'{e{Vh-V)~h), (1.5b) 

m' ^ {cTrn{V)y\e{V ^Vra)~m), (1.5c) 

where the boundary conditions are given by 

V{-oo)=Va, h{-oo) = l, m(-cx)) = 0, (1.6a) 

V{+oo) = V^, h{+oo) = 0, m{+oo) = 1. (1.6b) 



Here Va and Vui arc the pre- and post-front voltages, and Va < Vh < Vm < V^i- 
Equations (1.5) represent a system of fourth order so its general solution 
depends on four arbitrary constants. Together with constants 14, 14 and c 
this makes seven constants to be determined from the six boundary conditions 
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Fig. 1.3. (a) Tlie AP potential and (b) the gating variables h and m as functions 
of the travelling wave coordinate Z = z\/T). The solution of the CRN model [11] is 
given by broken lines and that of the three- variable model of (1.5) by solid lines. The 
prefront voltage and the excitation parameter in (1.5) are chosen as Va = —81.18 mV 
and j = 0.956, respectively and correspond to the equilibrium values in the realistic 
model. The gates h and m are indicated in the plot. The iterative analytical solution 
(1.14) is indistinguishable from the numerical solution of (1.5) after 30 iterations. 



in (1.6). Thus, we should have a one-parameter family of solutions, i.e. one of 
the parameters {Va, c) can be chosen arbitrary from a certain range. 
Comparison of the three- variable model (1.5) with the realistic CRN 
model [11]. The simplified three-variable model (1.3) and its ODE version 
(1.5) provide an excellent approximation to the fronts of the action potential 
in human atrial tissue as demonstrated in figure 1.3 where a comparison with 
the solution of the realistic CRN model [11] is presented. As must be expected 
the voltage in the simplified model remains constant after reaching its post- 
front value Vuj while the voltage in the CRN model assumes a shape typical 
for an action potential. To quantify further the comparison between the two 
models below we list the values of the wave speed, the post-front voltage 
and the maximum rate of AP rise. For the realistic model [11] these values 
are C = 0.2824 mm/ms, = 3.60 mV and idV/dt)„,ax = 173.83 V/s. 
The respective values for the simplified model (1.5) are C = 0.2372 mm/ms, 
Vcj = 2.89 mV and {dV/ dt),narc = 193.66 V/s. The relative errors made by 
the simplified model in estimating the wave speed and the the maximum rate 
of AP rise are 16% and 11%, respectively, and the absolute error in estimating 
the post-front voltage is —0.7 mV. 

We recall that our main motivation for the derivation of the three- variable 
simplified model was to reproduce the realistic front dissipation behaviour 
of atrial tissue as it appears for example in the CRN model [11]. The third 
column of figure 1.1 illustrates the dissipation of a propagating front in equa- 
tions (1.3) in response to a temporary block of excitability. It is observed that 
the dissipation behaviour of the simplified system resembles the one of the 
realistic CRN model In this aspect of the behaviour our model is superior to 
the simplified models of FitzHugh-Nagumo type. 
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1.3 Iterative analytical solution 

Exact solution for V < V„i. For V < V„x the reduced model (1.5), (1.6a) 
has a two-parameter family of exact solutions, with parameters and c. 
Since the boundary condition m(— oo) = is an equilibrium point of (1.5c), 
m(z) = remains a solution for all z < ^. It follows that (1.5a) is a constant- 
coefficient linear homogeneous equation in this interval, and its solution 

y{z)^V^ + {y,,^V^)e-\ (1.7) 

satisfies boundary conditions V^(— oo) = Vq, T^(0) = Vh and V^(^) — Vm, 
provided that the internal boundary point ^ is located at 

C-lnf^V (1.8) 



Finally, h{V) = 1 is a solution of (1.5b) with boundary condition h{Va) = 1 in 
the interval V < Vh, because it belongs to its equihbrium set. In the interval 
y > Vh, (l-5b) can be re-written in the form 

-^{\nh) = - \ (1.9) 

and its solution can be immediately obtained, 

dV 



W^exp(--1/J 



(1.10) 



Approximate solutions for V > Vm- For V > Vm, we rewrite (1.5) as 



dz \ dz 

-^{lnh) = ^— (1.11b) 

dz CTh(Z) 

where 

f{z) = -7^,{V{z)) J h{z) m3(z). (1.12) 
The boundary conditions are 

V{i) = Vn, (1.13a) 

V'{i)=c{Vm-Vo:), (1.13b) 

h{£,) = ho = exp ( -1 f I , (1.13c) 
V Jvu {V-Va)Th{V)l 
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m(e) = 0, (1.13d) 
V'{oo)=0. (1.13e) 



This problem is equivalent to the following system of integral equations 




The last equation Eq. 1.14d imposes an additional relationship between pa- 
rameters Va and c, so the ultimate solution depends only on one arbitrary 
parameter, say Va, which in reality may be determined by the pre-history of 
the medium through which the excitation front propagates. 

The system Eq. 1.6 can be solved by iterations, starting from a suitable 
initial approximation. The iterations converge for a number of various initial 
approximations, see figure 1.4. In section 1.4 we discuss two selected initial 
approximations. The first of them, Al, is favourable from a numerical point 
of view. The second one, A2, is important in the context of our recent work 
[30] since it leads to an even further formal simplification of problem (1.5) to 
a system which allows exact solution and extensive analytical study. 



1.4 Selected initial approximations 

Al: The small-diffusion initial approximation. A simple initial approx- 
imation may be obtained by considering a space-clampcd version of (1.5) cor- 
responding to the limit D — > of very small constant of diffusion in equations 
(1.2) (note that this is applied only to the solution in the interval V > Vm) and 
replacing the gating variable m with its quasi-stationary value m which in our 
asymptotic limit and for V > Vm equals 1. Hence, an initial approximation 
may be chosen to satisfy 

^=hr.{V)jh, (1.15a) 
^ = -h/{r,{V)). (1.15b) 
This system can be solved in quadratures. 
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Fig. 1.4. Convergence of the iterative solution (1.14), starting from initial approx- 
imation Al (triangles), same as Al but with the equation for m gate retained 
(squares) and A2 with V^"^ = Vm (circles). The excitation parameter and the pre- 
front voltage are j = 0.9775 and Va = —81.18 mV. The values on the right-hand 
side y-axis represent the numerical solution of the boundary value problem (1.5). 



dV 

KV)^ho~ (1.16a) 
/ , (1.16b) 

jv, iN.{v)jhivy 

where the initial conditions Vb and ho are given by Eq. 1.13a and Eq. 1.13c. 
A2: The 'caricature' initial approximation. An even simpler initial ap- 
proximation is 

1/ = const. (1.17) 

In this case functions of voltage /Na(l^), Th{V) and Tm{V) take constant values 
^(1^^°^), T-/,(y{o>) and T„(y{°>), respectively. Then quadratm-es (1.14) for 
the results of the first iteration are obtained in explicit formulae; moreover, 
(1.14d) is resolved explicitly. These explicit formulae have been reported in 
our recent work [30] (see formulae (9) of that paper). There this piecewise- 
linear simplification was considered only as an arbitrary "caricature" with 
the purpose of merely analysing qualitative features of the solution set. Here 
we note that it actually appears naturally as a step the iterative procedure 
leading to a numerically accurate solution. 

The iterative analytical solution (1.14) obtained from these initial approx- 
imations, is indistinguishable from the numerical solution of (1.5) after some 
30 iterations and has the shape of a travelling front as shown in figure 1.3. 
Convergence and uniqueness of the iterative solution (1.14). The 
iteration procedure produced by (1-14) is nonlinear and non-monotonic, and 
we do not have a rigorous proof of its convergence from any given initial ap- 
proximation. Likewise, we do not have a rigorous proof that the solution of 
the boundary-value problem (1.13) is unique. However it is straightforward 
to see that if the iterations converge, the result is a solution of the boundary 
value problem, due to the above mentioned equivalence of (1-13) and the sys- 
tem of integral equations (1-14). Figure 1.4 illustrates that the first several 
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iterations obtained from different initial approximations oscillate about the 
correct numerical solution of the problem and ultimately converge to it. The 
"small diffusion" initial approximation Al is particularly interesting because 
the results of the second iteration are already very close to the accurate nu- 
merical values for the wave speed and the post-front voltage. Indeed, if the 
excitability parameter and the pre-front voltage arc chosen at their physio- 
logical resting values j = 0.9775 and = —81.18 mV, respectively [11], at 
the second step of the iterations the value of the dimensional wave speed is 
C = 0.2255 mm/ms which has only 20% relative error compared to the value 
0.2824 mm/ms of the realistic ionic model (1.2) and 5% relative error com- 
pared to the value 0.2372 mm/ms obtained numerically from equations (1.5). 
Similarly the post-front voltage = —7.22 mV compares well with the value 
3.36 mV of the CRN model (1.2). Clearly, the second iteration obtained from 
initial approximation Al introduces certain errors. Numerically, however, it 
is immensely superior to any other numerical scheme because it involves only 
a single evaluation of formulae (1.16) and a two-fold evaluation of formulae 
(1.14). In addition, the dangers of numerical divergence associated with many 
of the alternative numerical schemes for solving problems (1.2) or (1.5) are 
avoided since the above expressions are mathematically well-behaved. 

The two initial approximations discussed above are essentially different. 
In the small-diffusion approximation the initial guess for the voltage V{z) 
is a function V'^ = v{z) while in the 'caricature' approximation it is a con- 
stant =const. The fact that two essentially different initial approxima- 
tions give the same limits is a strong indication that the iterative proce- 
dure leads to a unique solution. To support this claim further we have per- 
formed calculations starting from initial approximation A2 with initial values 
V^^^ ~ —28, —30, —32, —34, —36, —38. In all cases the iterations converged to 
the same solution, in a similar manner to those shown on figure 1.4. The 
small-diffusion initial approximation is not discussed here because its initial 
guess V'^ — v{z) is the unique solution (1.16) of equations (1.15) and thus it 
does not depend on any arbitrary parameters. 

So, although we cannot exclude the possibility that the iteration proce- 
dure may not converge from some "bad" initial guess, the examples considered 
provide an evidence that a reasonable initial approximation always gives con- 
verging iterations, and the solution is unique. 

1.5 Discussion 

We have presented an analytical approach to the description of the speed and 
the structure of an excitation front in a model of human atrial tissue [11]. We 
have identified small parameters in the realistic model and used asymptotic 
arguments to obtain a simplified three- variable model of the excitation front. 
Although we have explicitly used certain quantitative features of atrial tissue 
[11], the main properties used are generic for cardiac excitation models so 
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the approach should be applicable, possibly with suitable modifications, to 
other cardiac equations models too. Our model takes the form of a nonlinear 
eigenvalue problem with a piece- wise right-hand side defined over three voltage 
intervals. This model is solved explicitly in the first two intervals, and in the 
last interval we have suggested an analytical iterative procedure capable of 
producing a solution with a good accuracy already after the first iteration. 
The iterative procedure can be started from reasonably chosen simple initial 
approximations and converges to a unique solution which differs only within 
few percent from the solution of the realistic ionic atrial model. 

An important feature of our approach is that it is capable of correctly 
describing the excitation propagation at reduced excitabilities, up to a com- 
plete block of propagation via "front dissipation" mechanism, which is com- 
pletely unachievable by traditional analytical approaches based on FitzHugh- 
Nagumo type of equations. This aspect has been analysed in our earlier pub- 
lications [5, 6, 9, 30]. In particular, our recent study [30] is devoted to a 
discussion of one practical application of our approach. There we have used 
a numerical solution of the simplified three-variable model (1.5) to propose a 
simple criterion for break-up and self-termination of spiralling waves and have 
confirmed our predictions by numerical simulations of the realistic model of 
Courtemanche et ai. However, the important question of finding an analytical 
solution of our simplified model (1.5) was now solved in this paper. 

The possibility of obtaining numerically reasonable analytical approxima- 
tions to front solutions in realistic cardiac equations, demonstrated in this 
paper, opens the way for analytical description and, possibly, a better un- 
derstanding, of more complicated regimes in excitable media, such as wave 
break-ups and spiral waves. 

This study is supported by EPSRC grant GR/S75314/01. 
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